Finite temperature ordering of dilute graphene antiferromagnets 
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We employ large-scale quantum Monte Carlo simulations to study the magnetic ordering transition 
among dilute magnetic moments randomly localized on the graphene honeycomb lattice, induced 
by long-ranged RKKY interactions at low charge carrier concentration. In this regime the effec- 
tive exchange interactions are ferromagnetic within each sublattice, and antiferromagnetic between 
opposite sublattices, with an overall cubic decay of the interaction strength with the separation 
between the moments. We verify explicitly, that this commensurability leads to antiferromagnetic 
order among the magnetic moments below a finite transition temperature in this two-dimensional 
system. Furthermore, the ordering temperature shows a crossover in its power-law scaling with the 
moments' dilution from a low- to a high-concentration regime. 



I. INTRODUCTION 

Since the discovery of the single-layer hexagonal car- 
bon sheets of graphene, its fascinating physical proper- 
ties have stimulated a lot of research on the static and 
transport properties of fcrmions on the two-dimensional 
hexagonal lattice [1, 2]. Soon after its discovery, the role 
of disorder effects on graphene has opened up new re- 
search horizons [3, 4]. In particular, disorder induced 
localized states [5] and magnetism [6-9] have attracted a 
lot of interest. For instance, it was found that graphene's 
electronic properties give rise to the efficient forma- 
tion of magnetic moments from adatoms located on the 
graphene surface [8], or simply from defects [6, 7, 9]. Mo- 
tivated by such studies, we here consider a situation such 
as shown in Fig. 1, where local moments are induced in a 
graphene sheet by for instance adatoms or defects, each 
single magnetic moment being associated with a partic- 
ular lattice site of the underlying lattice [10, 11]. At low 
carrier concentration, the low density of states near the 
Dirac points suppresses the Kondo effect and allows for 
the formation of magnetically ordered states by the in- 
teraction between the localized moments and the conduc- 
tion electrons [10-12]. Namely, they induce long-ranged 
Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange in- 
teractions between the localized magnetic moments, de- 
scribed by the Hamiltonian 



RKKY 
ij 



(1) 



The effective interaction J^^^^ , mediated by itinerant 
electrons, strongly depends on the electronic properties 
at the Fermi energy. While in typical metals an oscil- 
lating coupling at 2k p is expected (kp being the Fermi 
wave vector) , decaying as 1 /r^ in two dimensions (r be- 
ing the relative distance between impurities), the semi- 
metallic properties of graphene lead to a different behav- 
ior [6, 10, 11]. Indeed, the absence of an extended Fermi 
surface leads to a cancellation of the 1/r^ term and leaves 
the next term decaying as l/r^ without any 2kp oscil- 
lations. Furthermore, it was revealed in Refs. [10, 11] 




FIG. 1: (Color online) Magnetic moments localized on a hon- 
eycomb lattice. Due to the commensurate nature of the long- 
ranged exchange interaction, the moments order antiferro- 
magnetically below a finite transition temperature. 



that pair-wise interactions are ferromagnetic (antiferro- 
magnetic) among the same (different) sublattice on the 
bipartite honeycomb lattice of graphene. 

This means for the exchange couplings in Eq. (1), that 



J. 



RKKY 



= e,jj(|r, -r^l). 



(2) 



with J > and e,- 



T (-f 1) if i and j belong to the 



same (different) sublattice [10, 11]. Here, r,; denotes the 
(random) position of the i-th magnetic moment Si on 
the honeycomb lattice, which we consider to be spin-i 
quantum spins, in order to explore the extreme quantum 
limit. Larger spin values of the moments will not qual- 
itatively change the results. The commensurate nature 
of the RKKY interactions was linked to the bipartiteness 
of the underlying lattice geometry [11]. In a more recent 
work [13], this general result was called into question for 
graphene nanoribbons, due to the presence of zero-mode 
contributions. In bulk graphene however, these correc- 
tions vanish in the thermodynamic limit [13], thus re- 
covering the commensurate form of the interactions in 
Eq. (2). In the following, we focus on the most basic 



2 



model that contains the main features of such exchange 
interactions (i.e. their commensurate, long-ranged na- 
ture) between the magnetic moments, and leave for dis- 
cussion at the end of the paper several aspects relevant 
to graphene, such as doping effects, a finite extension of 
the localized moments, and structural defects. 

The remainder of this paper is organized as follows: 
In the following section, we review some general results 
concerning long-range order in low-dimensional systems. 
Then, wc present our numerical approach in Sec. III. The 
results of our simulations are discussed in Sec. IV and 
Sec. V. Concluding remarks are made in Sec. VI, while 
an appendix provides details about the relevant length 
and energy scales on the diluted honeycomb lattice. In 
the appendix, wc furthermore introduce the notion of an 
effective coordination number for diluted magnetic mo- 
ments, that will be convenient for the discussion in Sec. 
V. 



II. LONG-RANGE ORDER IN 2D 

Starting from the effective exchange interactions of 
Eq. (2) in Eq. (1), wc explore its consequences for the 
finite-temperature ordering transition between magnetic 
moments on the honeycomb lattice. The stability of long- 
range magnetic order in d < 2 systems with power-law 
decaying interactions has been the subject of a large num- 
ber of theoretical studies in the past [14-27] . Regarding 
the Heisenberg model with 1/r" interactions, the semi- 
nal paper of Mermin and Wagner [14], proving the ab- 
sence of finite-T spontaneous order if a > d + 2, was 
recently reconsidered by Bruno [23], who gave stronger 
conditions, notably on the appearance of fcrromagnetism 
for oscillatory interactions. For instance, an interaction 
of the form cos(fcor)/r" (fco 7^ 0) cannot lead to finite- 
temperature fcrromagnetism if a > 5/2. For the case 
under study here, Bruno's result implies the absence of 
finite-temperature antiferromagnetic order for a > 2d 
which does not deviate from Mermin- Wagner's results in 
d = 2. 

Early renormalization group calculations on classical 
0(n) models with l/r'^+'^ couplings predicted [15, 16] a 
(T-dependent criticality with a finite ordering tempera- 
ture for a < d: Tc oc (d — <^)/in- — 1) when a ^ d [16]. 
For instance, while for a < d/2 exponents take exact 
"classical" values (77 = 2 — cr, ly = 1/ct, 7 = 1), our 
case (7 = 1 lies on the boundary of this regime where 
logarithmic corrections are expected [15]. In particular, 
the correlation length exponent turns out to be = l"'', 
which fulfills the Harris criterion [28, 29] v > 2/d. From 
such a statement we expect on general grounds that clean 
and disordered systems with power-law interactions dis- 
play similar critical behaviors if cr < 1. We indeed find, 
that the model in Eq. (1) exhibits at finite temperature a 
phase transition to an antiferromagnctically ordered Ncel 
state, both in case of a fully covered lattice and also for 
the case of diluted magnetic moments, with apparently 
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FIG. 2: (Color online) Temperature dependence of the stag- 
gered magnetization ruAF for different system sizes L = 
8, 16, 32, 40, 48, 56, and 64 for the occupation density p — 1 
(full coverage). The inset shows the temperature dependence 
of the specific heat Cv for system sizes L — 16, 48, and 64. 
The dashed lines indicate the position of the transition tem- 
perature as estimated from a finite size scaling analysis of the 
Binder parameter. 



similar critical exponents. 

In the following, we analyze in detail the dependence 
of the ordering temperature on the concentration p of 
the magnetic moments. While in the realistic parameter 
regime the magnetic moments will be dilute, i.e. p ^ 1, 
wc consider for completeness the whole range up to (and 
including) the case of full coverage p = I, where a mag- 
netic moment resides on every lattice site. For the full 
coverage case, we also performed simulations for an un- 
derlying square lattice, in order to explore more generally 
the magnetic ordering transition in quantum antiferro- 
magnets with long-range interactions in two dimensions. 
While in previous works [20-22], the case of solely fer- 
romagnetic interactions on a square lattice geometry has 
been analyzed, the current case and the efl:ects of dilution 
have not been considered thus far. 



III. METHODS 

For our study, we employed large-scale quantum Monte 
Carlo simulations based on the stochastic series expan- 
sion representation [30], with an improved diagonal up- 
date scheme adapted to systems with long-ranged inter- 
actions [31]. In addition, we used Walker's method of 
alias [32] in order to speed up the algorithm [33] . We per- 
formed simulations on finite systems with TV; = 2L^ lat- 
tice sites and linear system sizes ranging up to L = 192, 
depending on the concentration p of the magnetic mo- 
ments. For p < 1, we performed statistical averages over 
independent realizations of the moments' distribution 
in a canonical ensemble, such that no sample-to-sample 
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FIG. 3: (Color online) Data collapse of the Binder parameters 
for the full case p = 1 in a finite size scaling analysis. Here, 
t — {T — Tisi) /Tm denotes the reduced temperature. The inset 
shows the Binder parameters Q for different system sizes in 
the vicinity of the consecutive crossing points. 



fluctuations in the total number of moments N = p x Ni 
result. Typically, we performed disorder averages over 
several thousand realizations, and verified that the cal- 
culated observables followed Gaussian distributions, as 
expected. We always employed period boundary condi- 
tions, and used the minimum image convention for the 
1/r'^-decaying exchange constants. For each choice of p 
and L, we measured the staggered magnetization, ob- 
tained using the standard operator 



1 ^ 



after performing the disorder averaging as 



rriAF = \/[{S\p)]av 



(3) 



(4) 



Here, — ±1, depending on the sublattice to which 
spin i belongs on the honeycomb lattice, (...) denotes the 
QMC statistical mechanics expectation value for each re- 
alization of disorder, and [■.■]av the final disorder averag- 
ing. We also calculated the Binder parameter 



Q = [(sj, 



]at./[(S^_p) 



2 

av ■ 



(5) 



We then used a finite size scaling analysis to extract the 
critical temperature and exponents from the finite size 
data of iriAF and Q. More details about the finite size 
scaling analysis are provided below. 



IV. FULL COVERAGE 

As a useful starting point in the absence of any dis- 
order, we consider firstly the full coverage case p = 1. 



In Fig. 2, we present the QMC data for the temperature 
dependence of the staggered magnetization for various 
system sizes. This data shows, that a finite temperature 
ordering transition takes place near T ^ 1.3 J. This is 
also evident from the behavior of the specific heat Cy, 
shown in the inset of Fig. 2. It also exhibits pronounced 
finite size effects, that need to be accounted for in order 
to extract the transition temperature. 

For this purpose, we calculated the Binder parameter 
Q inside the transition region, for which the finite size 
data is shown in the inset of Fig. 3. The strongly moving 
crossing points in the data for consecutive system sizes 
again indicate significant finite size effects, that need to 
be taken into account in the further analysis. We thus 
employed a finite size scaling analysis including leading 
corrections to scaling [34] , used in several precision stud- 
ies on critical properties in quantum spin systems [35-37]. 

It is based on a scaling ansatz for the Binder parameter 



Q{t, i) = (1 + cL-^)gitL^^-' + dL-^'"), 



(6) 



with the reduced temperature t = {T — Tn)/T^, and 
the scaling function g. From this analysis, the critical 
exponent v is also obtained. Furthermore, c, and 
d describe the leading corrections to scaling, which are 
necessary in order to fit the QMC data obtained here 
for a system with long-ranged interactions on the limited 
system sizes available to our numerical study. Follow- 
ing Ref. 35, we represent g up to forth-order in a Tay- 
lor expansion {g{x) = go + gix + g2X^ + g^x^ + gix'^), 
and use bootstrapping in combination with a standard 
Levenberg-Marquardt nonlinear optimization algorithm 
to perform the minimization procedure and fit the nu- 
merical data to the above scaling from. 

In Fig. 3, we show the resulting data collapse for the 
case of p = 1. We find that the QMC data can be fitted 
well to the above scaling form, leading to an estimate 
of the transition temperature of T^/J = 1.298 ± 0.001 
with three significant digits. While the other fitting pa- 
rameters are less constrained by the finite-size data (see 
below) - as observed also in the above-mentioned high- 
precision studies of short-range interacting quantum spin 
systems [35-37] - we obtain from the finite-size analy- 
sis robust estimates of Tn , which is the quantity we are 
mainly interested in for this study; in particular, since we 
will analyse its dependence on the dilution p in the follow- 
ing section. Furthermore, we obtain an estimate for the 
correlation length critical exponent = 1.04±0.02. This 
value is in good agreement with the predicted value v — 1 
from the renormalization group approach [15], the small 
deviations from this prediction being attributed to loga- 
rithmic corrections in the dependence of the correlation 
length on the reduced temperature [15]. However, given 
the restricted range of system sizes available to our QMC 
study, we are not in a position, to accurately account for 
these additional corrections. From the same finite-size 
data, we also obtain estimates for w — 0.4 ± 0.1 and 
(f) = 0.6 ±0.1, which are less constrained within the boot- 
strapping analysis. Similar as for we expect residual 
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FIG. 4: (Color online) Temperature dependence of the stag- 
gered magnetization niAF for different system sizes L = 
32, 64, 96, 128, 160, and 192 for the occupation density p = 0.1. 
The dashed line indicates the position of the transition tem- 
perature as estimated from a finite size scaling analysis of the 
Binder parameters. 



FIG. 5: (Color online) Data collapse of the Binder parameters 
for p = 0.1 in a finite size scaling analysis. Here, t — {T — 
Tm)/Tm denotes the reduced temperature. The inset shows 
the Binder parameter Q for different system sizes taken in the 
vicinity of the consecutive crossing points. 



finite-size effects also on these values due to the logarith- 
mic corrections. These values are about a factor of two 
smaller than the values given e.g. in Rcf. 35 for the case of 
the quantum phase transition in bilayer Heisenberg mod- 
els, where however the expected universality class is that 
of the three-dimensional Heisenberg transition instead of 
the mean-field behavior expected here. From the fitting 
procedure, we obtain non-zero values for both prefactors 
of the subleading finite-size corrections, d = — 2±0.6, and 
c ~ — 0.11±0.06. This exhibits the necessity of including 
the subleading finite-size corrections to the leading scal- 
ing behavior in Eq. (6). The Taylor expansion coefficients 
of the scaling function g can be estimated from the boot- 
strapping analyis as go — 2.93 ± 0.04, gi — 0.18 ± 0.04, 
32 = -O.liO.Ol, 53 = -0.001±0.008, 34 = 0.009±0.003. 
The last two coefficients remain more unconstrained than 
the other fitting parameters. This indicates that g could 
also be represented well by a second order polynomial 
within the considered region close to T/v with coefficients 
similar to those given above. 

We performed an analysis for p = 1 also for the case of 
an underlying square lattice, and obtained the transition 
temperature in that case to be Tn / J = 1.855 ± 0.02, 
which is in fact close to the values obtained from 
QMC simulations and from using a Green's function 
decoupling for the fully ferromagnetic case [21, 22]. 
Furthermore, for the square lattice, we obtain an 
estimate of = 1.13 ± 0.07, which within the error 
bars agrees with the result for the honeycomb lattice, 
but deviates more from the mean-field value. From the 
finite size scaling of niAF oc L~^/'^ at T/v, we extract 
the ratio fi/v = 0.52 ±0.04, which within error bars 
agrees with the value 13/v = 1/2 from renormalization 
group calculations [15]. The estimates for lo = 0.4 ± 0.1 



and 4> = 0.6 ± 0.1, that we obtain for the square lattice, 
also agree within the error bars with the result for 
the honeycomb lattice. Again, this is expected, as the 
ordering transitions on both lattices belong to the same 
universality class. 



V. RANDOMLY DILUTED MOMENTS 

After having considered the full coverage limit, we now 
turn to the case of diluted magnetic moments, p < 1. 
Also in this case we do obtain a finite temperature order- 
ing transition. For example, the QMC data for the stag- 
gered magnetization at p = 0.1 is shown in Fig. 4. The 
corresponding data for the Binder parameter is shown in 
the inset of Fig. 5. Performing the same finite size scaling 
analysis as before, we estimate the Ncel temperature as 
Tn/J = 0.0357 ± 0.0006 for p 0.1. The correspond- 
ing data collapse of the Binder parameter is shown in 
the main panel of Fig. 5. The estimate for the cor- 
relation length critial exponent i' = 1.0 ± 0.02 appears 
somewhat closer to the mean-field value, while the results 
for w = 0.4 ± 0.1 and = 0.6 ± 0.1 agree well with the 
above values at p = 1. 

Proceeding in the same way for various values of p, 
we eventually obtain the dilution dependence of the Neel 
temperature shown in Fig. 6, which summarizes the main 
results from our numerical study. Concerning the esti- 
mates for the exponents u, w and 4>, we cannot observe, 
within statistical errors, any systematic changes with p 
from their values in the clean limit, which appears con- 
sistent with the discussion in Sec. II. 

On the other hand, the transition temperature shows 
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FIG. 6: (Color online) Magnetic transition temperature Tm 
as a function of the occupation density p. The inset shows 
the same data on a log-log plot. A power-law dependence 
proportional to p^^^ for p < 0.2 is indicated by the dashed 
line. 



also with the help of a p-dcpendcd effective coordination 
number Zesi^p), as defined in Eq. (A4), which displays 
two distinct regimes: For large dilution p <C 1 (i.e. for 
(r) 3> 1), ZcS ^ 1, and the natural energy scale for the 
magnetic ordering is set by the coupling value at the 
average distance (i.e. Jtyp), because each moment has 
only a few neighbors to couple with and thus Tn ^ Jtyp- 
Increasing p, once Zcs becomes significantly larger than 
one (i.e. once (r) ^ 1) the relevant energy scale which 
controls the ordering of the moments will be controlled 
by the average J^g, directly proportional to p. As 
shown in the appendix, the crossover between these 
two regimes takes place near p ^ 0.25. It is interesting 
to compare such a concentration to the percolation 
threshold of the 2D honeycomb lattice p* ^ 0.697 [38] 
where nearest-neighbor interacting quantum spins lose 
long-range magnetic order in the ground state [39]. The 
absence of any feature at p* in the present study is in 
fact consistent with the sizeable value of the effective 
coordination number Zcs{p*) 6.5. 



a strong dependence on p that we now analyse further. 
In the range 0.3 < p < 1, this dependence is almost per- 
fectly linear. An extrapolation of the linear suppression 
would exclude finite-temperature magnetic order below 
p 0.18. However, we find the low-p behavior of T/v to 
deviate from this linear behavior. In fact, as seen from 
the inset of Fig. 6, which shows the same data on a log- 
log plot, below p ^ 0.2, T/v exhibits an algebraic increase 
with p, scaling as 

Tjvcxp^/^ p<0.2. (7) 

In the following, we discuss the relevant energy scales 
behind the different behavior of Tn at high and low con- 
centration of the magnetic moments: 

On a two-dimensional lattice, dilute randomly dis- 
tributed magnetic moments are separated by an average 
distance that scales as (r) oc p~^^^ (on the honeycomb 
lattice (r) = p^^^'^/2, cf. the appendix) which defines a 
typical coupling strength Jtyp = J{{r)) oc p'^/^. In the 
low-p regime, we thus find that the Neel temperature 
scales with the characteristic energy scale set by Jtyp- 
At higher concentrations, the scaling of T/v with p be- 
comes more mean-field-like, namely directly proportional 
to the mean-field average coupling (see Eq. (A3) in the 
appendix) J^g ex p. This leads to the linear behavior in 
T/v observed at higher values of p. In this regime, the 
average nearest-neighbor distance between the magnetic 
moments is (r) ^ 1, and does not vary much as a func- 
tion of p. Its main effect is the reduction of exchange 
paths, as the number of bonds that each moment is asso- 
ciated with reduces linearly with p. The crossover results 
near p 0.25, corresponding to a concentration regime 
beyond which the average distance becomes (r) ~ 1, as 
shown in the appendix. 

The behavior of T/v can be qualitatively understood 



VI. DISCUSSION AND CONCLUSIONS 

Motivated by recent results on the properties of RKKY 
interactions between localized magnetic moments on 
graphcne [10, 11], we performed a systematic study of the 
finite-temperature ordering transition of dilute spin-1/2 
magnetic moments on the honeycomb lattice, induced by 
a commensurate long-ranged exchange interaction. We 
found that in the low dilution regime, where the effective 
coordination number is close to unity (i.e. the average 
separation (r) ^ 1), the Neel temperature scales with the 
typical coupling Jtyp (xp^^^. For larger occupations, the 
behavior crosses over to a mean-ficld-like linear reduction 
of the Neel temperature from its value in the full coverage 
case. We also presented estimates for the critical expo- 
nents /3 and i^, which within statistical errors are consis- 
tent with the prediction from previous renormalization 
group calculations for the ferromagnetic classical 0(3) 
model, given that additional logarithmic corrections are 
expected [15]. In our analysis, we considered the extreme 
quantum limit of S' = 1/2 magnetic moments. However, 
the physical picture will not change except that the Neel 
temperature will scale with S{S + 1) for higher quantum 
spins. For the future, it will be interesting to explore 
the critical properties of such diluted quantum magnets 
with long-ranged exchange interactions in more detail, 
also considering other decay rates of the exchange inter- 
actions. This would require the consideration of signifi- 
cantly larger lattices. Our main focus here was on the di- 
luted case, relevant to the physical situation in graphene. 

In the case of graphene, the RKKY coupling, con- 
trolled by the ratio between Coulomb repulsion U and 
band-width M^, is J ~ W^/W - 1 eV [6] which for a 
moderate concentration p ~ 10~^ would give a critical 
temperature ^ lOK. Of course this estimate is based on 
a very simple model of localized point-like magnetic im- 
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purities. A more realistic description should be able to 
incorporate (i) the spatial extension ^ of the defects, (ii) 
the holes/electrons doping effects, (iii) lattice distortions 
(ripples for instance). Regarding (i), a finite area for 
a defect is expected to move the crossover concentration 
p ~ 0.25 above which MF behavior Tn ^ p occurs to- 
wards a lower value p ^ 0.25/^^. (ii) As already discussed 
in Rcf. [6], holes/electrons doping shifts the Fermi energy, 
thus leading to a finite Fermi wave vector kp ^ ^frTc (ric 
being the carriers concentration). RKKY interactions 
will oscillate with a wave length \p ^ while 
the average distance between moments is (r) ~ 1/ ,Jp. 
Therefore the above analysis, which ignores 2fci? oscil- 
lating terms is expected to be valid provided ric ^ p. 
Alternatively, one expects the Neel order to be destroyed 
upon carrier doping in graphene sheets. For instance, 
using an electric field to control the Fermi level would 
render it possible to induce a transition from the Neel or- 
dered regime for kp ^ ^/p onto a more complex regime 
at fci? ^ yjp where competing interactions, i.e. mag- 
netic frustration, associated with random dilution are 
expected to provide all the ingredients to achieve spin- 
glass physics. We note that in the commensurate case 
at half-filling, the ferromagnetic and the antiferromag- 
netic exchange interactions actually have different pref- 
actors [fl]. However, this does not lead to any frus- 
tration, and hence including these prefactors will not de- 
stroy the finite temperature antiferromagnetic state, (iii) 
With respect to lattice distortions, it would be interesting 
to account for the characteristic ripples in the graphene 
structure [1 , 2] and explore its consequences on the mag- 
netic order, given the long-ranged nature of the exchange 
interactions. This would extend a recent study that con- 
sidered this interplay between structural and magnetic 
properties within an effective Ising model with exponen- 
tially suppressed exchange interactions on the order of 
several lattice spacings [40]. 



Two directions appear feasible to experimentally probe 
for the two-dimensional magnetism in graphene at finite 
temperatures: using magnetic adatoms like Mn for in- 
stance, or extrinsic defects [42] that could be created 
by irradiation. In addition to randomly distributed mo- 
ments, it will be interesting to explore the situation con- 
sidered in Ref. 8, where the magnetic moments are placed 
using STM techniques onto specific lattice sites, and to 
examine the magnetic states induced by the RKKY in- 
teractions. For such studies, the effects of frustration 
could lead to exotic magnetic phases, the study of which 
is however beyond the scope of the quantum Monte Carlo 
approach, due to the infamous sign-problem [41] . In that 
respect, future experiments on graphene might even be 
employed as a quantum simulator for such magnetic clus- 
ters. 
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Appendix A: Energy scales on the diluted 
honeycomb lattice 

In order to gain insight into the role played by various 
energy scales, we performed a numerical analysis on a 
L ~ 5000 diluted system, introducing a fraction p of 
magnetic moments randomly on the honeycomb lattice. 
The nearest-neighbor distance (i.e. the distance from a 
randomly chosen moment to the closest other moment) 
obeys a probability distribution (see the inset of Fig. 7), 
that at low concentrations p is very well described by 

P{r) ~ 2'Krp cxp{~Trpr^) , (Al) 

thus resulting in an average nearest-neighbor distance 
(r) = p~^/^/2, shown in Fig. 8. This leads to a typical 
coupling strength Jtyp = Jii^)) oc p^^'^. It is interesting 
to compare this to the average nearest-neighbor coupling 
J"" . defined as 

avg ' 

JZ, = / Jir)P{r)dr, (A2) 

which, at low concentration p ^ I, turns out to be (i) 
much larger than Jtyp and (ii) a linear function of p. On 
the other hand, the mean-field average coupling 

'^a'^g ^^E^d'-'-r^-l) (A3) 

compares well to J^g at low doping. But while J^g 
remains linear (J^g = 27rahGxC(3)p [43]) as p increases, 
J"^g approaches Jtyp for larger p. This p-dependence of 
the different energy scales is shown in Fig. 7. 
The effective coordination number, defined as 

ZaS = J^l/ Javg: (A4) 

clearly traces these two different regimes. For (r) ^ 1 
(i.e. beyond p ^ 1), the system is highly diluted and 
ZcS ^ 1 increases only slightly with p, whereas once (r) ~ 
1 (i.e. beyond p ^ 0.25), the effective number of magnetic 
neighbors increases much more rapidly, proportional to 
p. This difference in behavior directly follows from Fig. 8. 
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FIG. 7: (Color online) Dependence of the various energy 
scales on the occupation density p. The inset shows the his- 
togram of the nearest-neighbor distance obtained from simu- 
lations of a 1/ = 5000 system, along with the analytic formula 
for P{r). 




FIG. 8: (Color online) Average distance (r) between nearest- 
neighbor magnetic moments, computed (green squares) on a 
L = 5000 honeycomb lattice with a concentration p of mag- 
netic moments. The continuum formula (2y^)^^ (black line) 
holds at low doping. In this regime, the effective coordina- 
tion number ZeS, computed over the same sample (blue cir- 
cles), remains close to 1 and increase faster only for larger p, 
with the mean-field behavior (red dashed line) reached beyond 
p~0.25. 
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